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We analyze data from simulations of two- and three-dimensional (2D and 3D) glass-forming liquids using 
a correlation function defined in terms of a memory function with a negative inverse power-law tail. The self- 
intermediate function and the autocorrelation functions of pressure and shear stress are analyzed; the obtained 
fits are very good, at least as good as with a stretched exponential. In contrast to the stretched exponential, 
the key shape parameter — the exponent of the power-law tail — seems to be the same for all three correlation 
functions. It decreases from a value around 2 at high temperature to a value close to 1.58 (2D), 1.50 (3D) 
at low temperatures. At the same time the amplitude of the tail increases towards towards a limiting value 
corresponding to a diverging relaxation time, which is related to anomalous diffusion. On the other hand, 
careful analysis of the long time behavior in the case of the intermediate scattering function suggests that the 
memory function is cut-off exponentially, which avoids the divergence of the relaxation time. Repeating the fits 
with an exponential cut-off included indicates that the power-law exponent is in fact independent of temperature 
and close to 1.58/1.50 over the whole range. Instead of the divergence, a fragile-to-strong crossover in the 
dynamics, estimated to occur around T = 0.40 for the 3D Kob-Andersen system. Another key parameter of 
the fitting procedure may be interpreted as a short-time rate, the amount of decorrelation that occurs in a fixed, 
relatively short time interval (compared to the alpha time). This quantity is observed to have a near-Arrhenius 
temperature dependence, while its wavenumber dependence seems to be diffusive {q~^) over a wider range than 
that of the relaxation itself, a further indication that this "bare relaxation rate" is simpler than the full dynamics. 



I. INTRODUCTION 

Time-dependent correlation functions are a standard tool 
for quantifying the dynamics of physical systems.' For the 
most basic analyses, knowledge of a single quantity — the 
relaxation time — and its dependence on external parameters 
such as temperature, is sufficient. But often more detailed 
information, in particular concerning the shape of the corre- 
lation function, is required. This might be, for example, be- 
cause competing theoretical explanations make different pre- 
dictions about the shape, for example. In the case of highly 
viscous glass-forming liquids^'^ some functions typically in- 
vestigated include the autocorrelation functions of energy, 
pressure, shear-stress, and dielectric moment. These are re- 
lated through the fluctuation-dissipation theorem' to the cor- 
responding dynamic response functions: frequency dependent 
specific-heat, bulk modulus, shear-modulus (or equivalently 
viscosity) and dielectric constant, respectively. The key fea- 
tures of viscous-liquid dynamics are (1) the non-Arrhenius 
temperature dependence of the main, or alpha, relaxation time 
and (2) the non-exponentiality of the corresponding coiTela- 
tion functions. Any theory of glass-forming liquids has to be 
able to explain these "nons" — violations of what is generally 
expected in relaxing systems (Arrhenius temperature depen- 
dence and exponential, or Debye, relaxation). 

This paper presents a method for fitting relaxation data in 
the time domain, by parameterizing not the correlation itself 
but its associated memory function."* The latter is a unique 
function related to the coiTelation function. The motivation 
comes from the hope that a relatively simple description of 
the memory function could be attained, and this hope turns 
out to be justified — very good fits can be achieved by taking 
the memory function to have the form of an inverse power law 
for non-zero times. Moreover the exponent of the power law 



seems to have a common value, approaching at low temper- 
atures a value around 1.6 for the 2D system and 1.5 for the 
3D system, for the different coiTelation functions, something 
that is not true of the "stretching parameter" (3 when fitting 
using a stretched exponential. In fact when the memory func- 
tion parameterization is generalized to include an exponential 
cut-off, the exponent is temperature independent. 

In the next section we discuss some general ideas for de- 
scribing relaxation and review the concept of the memory 
function associated with a correlation function. Following 
that is a brief description of our simulations, while SectionlTVl 
presents the details of our analysis method as applied to time- 
domain data from the simulations. Results from the fitting 
procedure are presented in Section [V] Section [Vl] discusses 
possible ways of interpreting the power-law description. 



II. NON-EXPONENTIAL RELAXATION AND THE 
MEMORY FUNCTION 

A. Correlation functions 

In this work we consider normalized coiTelation functions 
4'{t), typically (but not exclusively) defined in terms of a dy- 
namical variable A{t), as 



{AA{0)AA{t)) 
((AA(0))2} ■ 



(1) 



For example A{t) could be the potential energy, pressure, or 
shear stress. Here A indicates deviation from the thermal av- 
erage. Debye relaxation corresponds to a simple exponential 
decay with time constant tu: 
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Ipoit) = CXp(-t/Tl3). 



(2) 



There are many ways to define a relaxation time Tq, which 
agrees with tjj in the case of Debye relaxation. For our pur- 
poses the most convenient is the time integral of ^{t): 



i}{t)dt 



(3) 



One advantage of this definition is that it uses all of the in- 
formation contained in •4j{t), unlike, for example, a definition 
based on the time at which '4>{t) attains a certain value. It is 
also convenient when working with Laplace transforms, since 
this is simply the the s = value of the Laplace transform of 
ilj{t). The function most commonly used to fit time-domain 
data is the stretched exponential or Kohlrausch- Watt-Williams 
(KWW)^^'^ function. 



^Kww{t) = exp {-{t/TKwwf'''^'^) 



(4) 



where tkww sets the time scale and (3kww is a number be- 
tween and 1 , known as the stretching parameter Exponen- 
tial behavior is recovered for (3 — 1. By integrating we find 
that Ta = T{1/(3kww)tkww/Pkww- Fits to Eq. dH are 
often very good, but on the other hand experiments^ suggest 
that the relaxation crosses over to exponential at the longest 
times (a decade or more beyond tkww), corresponding to the 
last few percent of relaxation. This has not been investigated 
much, if at all, in simulations, probably because it is difficult 
to get high quality data at the longest times. In the case of 
single-particle diffusion, however, the crossover to exponen- 
tial behavior should coincide with the onset of Fickian diffu- 
sion and the convergence of the van Hove correlation function 
to a Gaussian, which have been studied by various authorsi^ 



B. The memory function 

The memory function concept was introduced by the work 
of Zwanzig'" and Mori," which provided a theoretical for- 
malism for the calculation of correlation functions of many- 
body systems, based on projector-operator techniques. One of 
the historically most important uses for the memory function 
has been to elucidate the short-time structure of correlation 
functions. For example, pure exponential decay, represented 
by a delta-like memory function, is unphysical at short times, 
since the correlation, being even and smooth, should have zero 
slope at i = 0. Replacing the delta function with a less sin- 
gular function (for example a step function, a Gaussian, or an 
exponential) gives more physical short-time behavior In this 
work, however, we are concerned with the long time behavior, 
and in our analysis will ignore the short time behavior, which 
is mainly associated with vibrational motion and irrelevant for 
structural relaxation. 

The memory function, K{t), is defined^ for a normalized 
autocorrelation function tp{t) by 



K{t - t')il){t')dt' . 



(5) 



The relation between tp{t) and K{t) is completely invertible, 
which can be seen by considering the Laplace transformation 
of this equation: 



sip{s) 



which implies 



K{sy 



(6) 



(7) 



an invertible relation between if^is) and K{s). Thus we can 
think of the memory function as a kind of transform of the 
correlation function. Setting s = in Eq. (|7]i gives an expres- 
sion for Tn '■ 



(8) 



This will be useful below. Although though of our anal- 
ysis will involve a pure power-law decay of K{t), it should 
be noted that if ijj{t) is known to decay exponentially at the 
longest times, the same must be true of the memory function: 
If ipit) decays as an exponentially ^ exp (—At) at long times 
(with A is real and positive), this implies that the least nega- 
tive singularity of ip{s) is at s — —A. In particular ip{s) is 
analytic, as well as real and positive, in a finite region about 
s — (recall ip{s = 0) = r^). The same is true of l/tp{s) 
and therefore K{s) {=l/ip{s) — s). Being finite and analytic 
at s = is incompatible with a power-law decay for K{t), 
which would correspond to a power-law for K{s) with (in 
general) non-integer exponent as s ^ 0. If an exponential 
cut-off is included in K{t), however, this moves the singular 
point a finite distance to the left along the negative real s-axis. 



C. Interpretation as a rate 

A useful starting point for considering the meaning of K{t) 
is the case K{t) — 'yS{t), which gives exponential relaxation 
i/j{t) — exp{—jt), i.e., Ta — 1/7, consistent with Eq. ([8j. 
This is the memory-less, or Markov, case. Note that the 'am- 
plitude' 7 gives the rate of the exponential decay. Consider 
integrating Eq. (|5]l as if it were a physical equation of motion: 
In the memoryless case, as time goes on, the delta function 
means that the only contribution to the integral on the right- 
hand side is from the current time, t' = t. Thus, the rate of 
change of ip{t) is proportional to its current value and inde- 
pendent of its previous values. This of course gives exponen- 
tial relaxation. If we generalize slightly to the case iir(t) — K, 
a constant, for t < T, and zero otherwise, then for times much 
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greater than T we again expect exponential decay with relax- 
ation time Tq, — 1/{KT)P When K{t) is non-zero for all 
t > 0, but /o K{t')dt' converges sufficiently quickly (e.g. 
exponentially) to a finite positive value, we can expect expo- 
nential decay in ip{t) for long times. 

For times t when the memory function is still changing, 
K{t')dt' gives an estimate of the instantaneous decay rate 
of 'ip{t). What happens when Kit), which must be positive 
at time zero (otherwise ^/'(i) will diverge), becomes negative? 
At that point Jp K{t')dt' starts to decrease, implying that the 
decay rate decreases. As long as K{t) < 0, ip{t) will exhibit 
slower-than-exponential relaxation — the instantaneous decay 
rate keeps decreasing. Therefore we expect that for viscous 
liquids the memory function is negative and significant for 
times up to Tq, around which time Jg K{t')dt' has more or 
less converged to the limiting value and exponential relaxation 
takes over. 



ifo = l-^i/V^o = l-V'i, (10) 

(since we assume the normalization i/jq — 1). Knowing now 
A'o, writing Eq. (|9]l for n = 1 gives again an equation with 
only one unknown, namely Ki . In this way the values of K,n 
can be solved for one by one. This process may be formalized 
using the idea of a generating function.'^ We consider the se- 
ries {V'ji} as the coefficients of a power series in a variable 
z, which (formally) defines the so-called generating function 
'^{z) = Y.'o V'mZ™, similai-ly K{z) = Y,"^ K^z"^- We may 
treat generating functions as quantities which can be added, 
multiplied and divided (providing the first coefficient is non- 
zero), and thus can note that the right hand of Eq. ^ is in fact 
the nth coefficient in the product series —K{z)4>{z) (this is 
like the convolution-product correspondence in Laplace trans- 
forms). Thus the equation may be re- written as 



D. Discrete memory function 

When we consider simulation data we will be not inter- 
ested in short-time vibrational contributions to ip{t) or K{t). 
A standard way of removing the effects of vibrations is to 
consider the so-called inherent dynamics,'^ obtained by mini- 
mizing configurations along the actual simulated trajectory to 
the corresponding local minimum of the potential energy — 
the inherent state. The correspondence is defined by steepest- 
descent minimization and is, in principle, unique for almost 
all configurations. It has been shown that the resulting tra- 
jectory of inherent states preserves the long time features of 
the dynamics.-i^^ Correlation functions are essentially un- 
changed except that the usual initial decay due to vibrations 
on the picosecond time scale is missing. A technical prob- 
lem is introduced by this procedure, however: The inherent 
correlation functions tend to have non-zero, and therefore dis- 
continuous, slopes at t = 0. This is not so surprising, given 
that the inherent trajectory is itself necessarily discontinuous, 
but it does pose problems for defining the short-time behavior 
of the memory function. To avoid these problems, and given 
that simulation data is generally available at regular, discrete 
times, we write a discrete version of Eq. (|5]l: 



V'n+l - "0" = ~ X! Kn-m1pm- (9) 
m=0 

Here time has been discretized in units of At, so 1/1,1 = 
= nAt). Equation ^ can be straightforwardly used to 
calculate the discrete values of the autocorrelation function, 
{i/n}, n > 0, given the discrete values of the memory func- 
tion, {Km], m > 0: Start by setting ?/;o = 1 and then for each 
71 > in turn, calculate the sum on the right hand side, which 
only involves "ipm for m < n and the known {Km}, to get 
the next unknown value, ipn+i- It is almost as straightforward 
to see that the reverse transformation is possible, i.e., given 
{i/n}, n > 0, to determine {Km}, m >Q. For example, from 
Eq. ^ for n — one finds 



{Aij){z) = -K{z)^/j{z), (11) 

where the difference function ( Ai/;) (z) is defined as the series 

(AV')o = 0, {Atp)n = ipn- ipn-i,n > 0. K (z) is then given 
by 

^W = -(^)(AV^)W: (12) 

which may be evaluated using standard algorithms for multi- 
plying and dividing power series.''' Mathematically this pro- 
cedure is equivalent to the inversion procedure mentioned 
above; thinking in terms of generating functions simply makes 
the implementation more straightforward. We use this trans- 
formation in the data analysis to relate correlation functions to 
their corresponding memory functions. In the following, we 
use K{t) to mean Km/ (Ai)^, where t — mAt. 

III. SIMULATIONS 

A. Systems and potentials 

We have simulated both two- and three-dimensional (2D 
and 3D) binary Lennard-Jones (BLJ) fluids. The parameters 
for the BLJ potential for both 2D and 3D are (where L and S 
stand for large and small particles, and e and a the energy and 
length scales, respectively) ell — 1, ess — 0.5, e^s — 1-5, 
(Tll ~ 1, fss = 0.88, (Tis = 0.8. All particles have the 
same mass m ~ 1. These parameters are identical to those of 
the 3D BLJ introduced by Kob and Andersen'^. The potential 
was truncated using an interpolating polynomial between 2.4 
aap and 2.7 {a, (3 G {L, 5}). In all simulations periodic 
boundary conditions with constant area/volume were used. In 
2D the total number of particles Np=100 of which 60% were 
of type L, while in 3D, iVp=1372 of which 80% were of type 
L. The particle density was 1 .2(7^^ in both cases (where d is 
the dimension). Constant energy dynamics were simulated us- 
ing the Verlet algorithm with a time step of 0.01 oj^i^^J mj tj^j^ 
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From now on, all quantities will be reported in the units de- 
fined by e^i, a^i^ and m. If one converts to physical units 
by assuming that the parameters for large particles are chosen 
to model Argon (i.e., aLL = 0.34 nm, ell — 997 kJ/mol, 
m = 39.95 u), then the time unit corresponds roughly to 2 ps. 



B. Runs 

Initial configurations consisted of fee lattices with parti- 
cles randomly assigned to be of type L or S such that the 
appropriate fraction was achieved. All simulations were at 
constant NVE; in particular the energy was controlled rather 
than the temperature T, but we quote the temperature defined 
in terms of the mean kinetic energy per particle instead as it 
is more meaningful. Starting at at T ^ 1.0 (where the liq- 
uid is not particularly viscous) ten independent runs (differing 
by the random placement of the particles on the initial lat- 
tice) were made. The initial lattices melted immediately and 
equilibrated in the liquid state. For each state point, an equi- 
libration run was carried out followed by a production run 
of the same length. Equilibration was checked by examin- 
ing the self-intermediate scattering function F{qm, t) for the 
larger particles; is 6.5 in 2D and 7.3 in 3D (corresponding 
roughly to the peak of the structure factor for large particles). 
F{qm,t) was required to be essentially the same for all ten 
runs, and the run-length to be of order 100 times the struc- 
tural relaxation time (r^ = /g F{q„i,t)dt). If the condition 
was not satisfied, the production run was considered to be part 
of the equilibration and a new production run was initiated 
from the final configurations. Steps in temperature were made 
by changing the total energy according to the estimated spe- 
cific heat; a new kinetic energy was chosen such that when 
summed with the current potential energy the correct total en- 
ergy would result. Then the particles' velocities were ran- 
domized using a Gaussian distribution with the appropriate 
width and then rescaling to give the exact desired kinetic en- 
ergy. During the equilibration process it was checked that the 
new temperature (mean kinetic energy) was the desired one 
and velocities were rescaled if necessary (generally a small 
change). 



C. Calculating correlation functions 

During the simulations, collective quantities such as poten- 
tial and kinetic energy and the different components of stress 
were computed and written at regular intervals (typically ev- 
ery 100 or 500 time steps, though some shorter runs were 
made with more frequent output in order to determine the 
short time behavior of the correlation functions). Configu- 
rations were saved at "logarithmic intervals" to allow calcu- 
lation of F{q,t) for a broad range of times t without writ- 
ing configurations as frequently as the collective quantities 
(which would require a huge amount of disk storage). Be- 
cause F{q, t) involves an average over particles, it has signifi- 
cantly smaller statistical error and therefore has low noise out 
to quite long times beyond the relaxation time. It should be 
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FIG. 1: (Color online) Left panel, shear stress (unnormalized) au- 
tocorrelation function for true and inherent dynamics. They differ 
only at short times. The slight dip in the curve for true dynamics at 
t = 5 is from combining separate short- and long- time simulations. 
Minimization was carried out every ten time steps, which defines the 
resolution of the memory function, At — lOdt = 0.1. The dashed 
line shows the correlation function obtained from attempting to in- 
vert the attempted fit to the memory function, while the dotted line 
shows that obtained by inverting using the exact memory function 
truncated alt = 100. Right panel, memory function for the inherent 
shear-stress autocorrelation function. Inset shows the negative tail in 
a double-logarithmic plot, along with a power-law fit. 



noted, however, that it is technically not the autocorrelation 
function of a dynamical variable. Autocorrelation functions 
of collective quantities were calculated using standard Fourier 
transform techniques. Both these and F{q,t) were averaged 
over the ten runs. For the main part of the analysis the auto- 
correlation functions were also "logarithmically averaged" in 
time. This weights different time scales equally and involves 
grouping the data points into equal-sized bins of log(i), and 
averaging within each bin. 



D. Inherent states 

Some of the analysis was caiTied on inherent-state con- 
figurations. The method used to minimize the energy was 
a combination of "molecular dynamics minimization"'^ and 
conjugate gradient."* This combination has been found to re- 
duce the possibility of finding the "wrong" minimum, which 
can happen if conjugate gradient alone is used.'^ Configura- 
tions were minimized every 10, 100 or 1000 configurations — 
more frequent minimization enables resolution of shorter time 
scales, but longer time scales cannot be accurately probed. 



IV. DATA ANALYSIS: DETERMINING K{t) 

A. Inversion of the correlation function to get the exact K{t) 

We start by applying the transformation described above 
to the simulation data. Fig. [T] shows the 2D, T=0.34 shear- 
stress autocorrelation function for both true and inherent dy- 
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namics. The agreement between the two at long times is clear, 
although there is more noise in the inherent case — due to the 
cost of minimization, not as long times can be simulated. The 
right panel of Fig. [T] shows the result of the discrete memory 
transform. Recall that this returns the exact memory function 
K{t) for the data. This is positive at t = 0, but negative for 
alH > (note that t really means nAt for an integer n; since 
minimization was carried out every ten time steps. At = 0.1). 
As shown in the inset, the data suggest an inverse power-law 
for the negative tail, although the noise is too great to justify 
this for t > 2, a very short time compared to the relaxation 
time Ta ^ 3000. Nevertheless let us attempt a fit to a simple 
power-law in the hope that this is in fact the correct form of 
K{t)fort > 0. The fit yields the expression 0.0512/iii86 for 

K{t). Testing whether this expression is in fact a good rep- 
resentation of the true memory function is simple: we simply 
use Eq. (|9|l to construct the correlation function correspond- 
ing to the fitted K{t). The value of Kq is not fit; rather the 
exact value from the inversion is used. The result, shown as 
the dashed line in the left panel of Fig. [T] is a disaster: af- 
ter reasonable agreement up to i ^ 10, the transformed ■i/;(t) 
stops decaying and in fact begins to increase and eventually 
diverges. The problem can actually be foreseen by consider- 
ing the time integral of K{t), which in the discrete case is the 
sum J2n At. The value of K[){— 1 — ipi) here is 0.0286, 
so that its contribution to the sum is 0.286. The sum over the 
power-law part requires the evaluation of the Riemann zeta 
function C,{s) at s = 1.186, which is estimated numerically to 
be 5.97. Including the appropriate factors of At and the pref- 
actor yields the contribution from the power law to be -0.469, 
making the whole sum negative. Moreover the sum first be- 
comes negative around where the inverted correlation function 
begins to increase-at this point the effective decay rate has 
become negative, and so the correlation function grows rather 
than decays. 

Clearly, direct fitting of the memory function is problem- 
atic. The fit which applies at least up to < ^ 2 clearly is 
wrong for t > 10. Because the noise in the exact memory 
function is so much larger than the true value, standard fitting 
procedures cannot distinguish it from zero in this region. On 
the other hand, we can directly check that it is indeed signif- 
icantly different from zero for quite long times: Taking the 
exact memory function, but setting it to zero for times greater 
than t = 100, and inverting, we find the correlation function 
indicated by the dotted line in Fig.[T] After t ^ 500 it drops 
below the true correlation function in a more or less exponen- 
tial fashion, as expected when the memory function is only 
non-zero for a finite range of time. Thus the apparent power- 
law behavior at short times must cross over to something more 
rapidly decaying, but still significantly different from zero. 



B. Matching the correlation function by fitting K{t). 

An approach other than direct fitting of the exact K{t) is 
clearly needed. Since the noise in the correlation function 
is small compared to the function itself, comparing this to a 
trial function in a fitting procedure will not suffer the same 




FIG. 2: Main figure, fitted memory function consisting of an inverse 
power-law multiplied by a piecewise linear function. Inset, compari- 
son of the shear stress autocorrelation function from simulations and 
the correlation function corresponding to the fitted memory function. 
In this case the correlation function was a combination of that deter- 
mined by inherent dynamics at short time and that determined by 
real dynamics at long-time, in order to have a more accurate and less 
noisy curve at long times. 

problem. This suggests a practical method for determining 
the memory function: We choose a functional form which is 
quite flexible, containing many parameters. For any given set 
of parameter values we can generate the corresponding trial 
correlation function which can be compared in a least-squares 
sense to the actual correlation function. Then we vary the pa- 
rameters to optimize the match. 

The value of Kq is fixed from the start to 1 — V'l- Given 
that the exact K{t) seems to have an initial inverse power- 
law decay for i > 0, we start with such a function. To pro- 
vide flexibility we multiply by a piecewise linear function; 
such a function will have discontinuous changes in slope but 
is simple and general — and hopefully the dominant behavior 
has been captured by the power-law factor. Thus this general 
fitting form Kg^n [t] is 

Kg,n{t) = -^L{t),t>0 (13) 

where L{t) is a piecewise linear function defined by a set of 
nodes tk and values Lk for fc = 0, 1, 2, . . . , Unodes- The first 
node is at t — 0, i.e., tQ — and Lq is fixed at unity. We 
choose Unodes ~ 20 to provide sufficient flexibility; a typical 
set of node locations tk is 0, 1, 2, 4, 5, 6, 8, 10, 20, 30, 40, 50, 
60, 80, 100, 200, 500, 750, 1000, 3000, 5000. The adjustable 
parameters are then the overall coefficient A, the exponent a, 
and the coefficients L^, k — 1, . . . , rinodes- At this stage we 
do not seek the true functional form of K{t); rather, we wish 
to have a fairly accurate representation which is relatively free 
of noise. To obtain the optimal parameter values in Kgen{t), 
a conjugate-gradient procedure is used. Some details of this 
are given in the appendix. Fig.|2]shows the result of this pro- 
cedure. While the discreteness due to the piecewise linear 
function is clearly visible, the overall form is clear. In particu- 
lar, there is a clear cross-over at around t = 5 from the power 
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FIG. 3: . (Color online) Comparison of autocorrelation function and 
fit with a memory function involving a single inverse power law in 
the tail. The data is the same as in Fig. |2l but only times which 
were multiples of 10 were included. The fitted exponent was 1.84, 
reasonably close to the apparent exponent of 1.7 observed in Fig.[2l 



law we initially identified to a more rapidly decaying power 
law, with exponent ^1.7. This cross-over means that the time 
integral converges to a positive value and the effective decay 
rate never becomes negative. 



C. Simple power-law form for K{t) for t > 10. 




t 



FIG. 4: (Color online) Main figure, examples of fits using the 3- 
parameter power-law form of K{t) to three correlation functions cal- 
culated from simulations of the Kob- Andersen binary Lennard- Jones 
mixture in 3D at r=0.46. Solid lines: data; dashed lines: fits. In con- 
trast to Fig. [3] logarithmic binning/averaging was used for the shear 
stress and pressure autocorrelation functions. 



Km = -Ko-^ — ,m>0, (14) 

where ({a) = X^j^il/"-" '^^e Riemann zeta-function. 
Then the relaxation time is given simply by 



The observation that K{t) is well approximated by a single 
inverse power law for times greater than about 5 immediately 
suggests a simplification: If we choose to not resolve times 
shorter than this, a fit using only a single power law for the 
non-zero part of K{t) may be possible. Fig.[3]shows that this 
is indeed true. Here times greater than 10 were included, cor- 
responding to sampling at an interval At = 10. This time 
scale is also significant for another reason (not necessarily un- 
related): It corresponds to the end of the vibrational contribu- 
tion to the stress autocorrelation function of the true dynamics 
(see left panel of Fig. [T]i. This suggests that we can avoid 
the time consuming energy minimization process altogether 
and use only the data from the true dynamics. There is one 
problem, though: We cannot get the t = value of the (in- 
herent) autocorrelation function from the true dynamics. To 
get around this we leave it as a fitting parameter; this gives 
us a procedure for fitting autocorrelation functions with three 
parameters: Rqi = 4>q/^i, A and a. Note that i?oi is directly 
related to Kq, so one could equivalently think of Kq as the pa- 
rameter. This fitting function involving a single inverse power 
law has been used in most of the analysis. 

It is convenient to represent that amplitude of the power- 
law in terms of a "tail-weight fraction" /, defined such that 
f — 1 corresponds to the limiting case where the sum over 
the whole negative (m > 0) part of A',„ exactly cancels the 
positive contribution from Kq- Thus we represent Km for 
TO > as 



^„(i^„/(Ai)2)Ai K„Atl-f 

This form makes explicit a few things: for near zero /, the 
relaxation rate is essentially the inverse of A'o- As / increases 
the negative tail "cancels out" part of the latter, increasing the 
relaxation time, and the limit / ^ 1 gives a diverging relax- 
ation time. 

Examples of fits to the three functions, F{q,n,t) and the 
shear stress and pressure autocorrelation functions, are shown 
in the main part of Fig. |4]for one temperature in the 3D sys- 
tem. Here and in the all the analysis presented below, the 
discretization interval was Ai=20. 



D. Long-time behavior: An exponential cut-off in K{t) 

As mentioned in Section|ll] experiments^ suggest that auto- 
correlation functions switch over to exponential decay at very 
long times. This is incompatible with the fitting function we 
have identified — as discussed in Section HH power-law decay 
of K{t) corresponds to a power-law singularity in K{s) as 
s ^ 0, implying non-analyticity of ■0(s) in the same limit, 
which excludes exponential behavior of K{t) at long times. 
Thus we should — in principle at least — allow the possibility 
of including an exponential cut-off in the fitting function for 
K{t). The question is whether the available data requires it to 
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FIG. 5: (Color online) Comparison of memory functions obtained for 
intermediate scattering function (3D, T — 0.46, Tq, — 1500) using 
3- and 4-parameter fits. The exponential cut-off in the latter is clearly 
visible. Inset, comparison of the correlation functions from the fit 
with the data, along with a stretched exponential. Here a double-log 
scale is used to emphasize the long-time behavior (the differences 
would not be visible otherwise). 
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FIG. 6: Arrhenius plot of relaxation times for structure (F{qm,t)), 
shear stress and pressure autocorrelation functions, for 2D and 3D 
data (closed symbol), determined by summing the fitted normalized 
memory function. 



V. RESULTS 



A. Relaxation times 



get a satisfactory fit. Only in the case of the intermediate scat- 
tering function is the noise sufficiently low to allow investiga- 
tion of the long time behavior (when the correlation function 
has decayed to less than 1% of its original value). Fig.|5]com- 
pares the three-parameter fit using Eq. [14] with a 4-parameter 
fit involving an exponential cut-off: 



^^^^,m>0 (16) 



where Tc is the characteristic time of the exponential cut-off. 
The memory functions are shown in the main part of the fig- 
ure, while the correlation functions are shown in the inset. For 
comparison, a stretched-exponential fit is also included. The 
three-parameter (pure power-law) memory-function fit notice- 
ably underestimates the decay rate at long times. Including 
the cut-off gives a much better fit, although some improve- 
ment is of course expected due to the extra parameter). Notice 
that the stretched exponential function fits also better than the 
pure power-law fit at long times, as well as the cut-off mem- 
ory function. For the shear stress and pressure autocorrelation 
functions the long time data is not good enough to warrant 
including the exponential cut-off as a fourth parameter, al- 
though we shall see that it can make sense to do so when there 
is reason to fix the exponent a to a particular value. In the 
following, "four-parameter" and "three-parameter" fits refer 
respectively to whether the exponential cut-off was included 
or not, regardless of whether one or more of those parameters 
was constrained in the fitting process. 



We start by presenting the relaxation times for the dif- 
ferent correlation functions. As stated earlier, we define 
as the integral of the normalized correlation function. With 
simulation data there is a contribution at short times due to 
vibrations. The relative height of this varies among the dif- 
ferent correlation functions, making the choice of normaliza- 
tion problematic. In principle this can be removed by con- 
sidering correlation functions of inherent quantities, but the 
time taken for quenches means that determining the correla- 
tion function accurately at long times is difficult. We avoid 
these problems by using the fitted correlation functions, nor- 
malizing by the zero-time value that emerges from the fitting 
process. To be consistent with the discrete-time formalism 
for dealing with the memory function we sum the correlation 
functions and multiply by At, rather than integrating them. 
This makes a positive difference of order At/2, negligible 
except when the relaxation time is of order At. This is the 
case at the highest temperatures, particularly for shear stress 
(note the upwards bend in the relaxation time curve at the ex- 
treme left of the plot). The value of Tq determined from the 
fit varies only slightly according to whether the 3- or 4- pa- 
rameter fitting function for K{t) is used. Fig.|6]shows Arrhe- 
nius plots for both 2D and 3D systems of the three coiTela- 
tion functions investigated. Some curvature — corresponding 
to so-called "fragility" of the viscous liquid — is evident, al- 
though less so in the 2D data. In both 2D and 3D the shear 
relaxation time is noticeably shorter than those of pressure 
and structural relaxation (F{q„i, t)). The pressure relaxation 
time tracks closely the structural one in 3D, but exceeds it 
noticeably in 2D. In fact the difference between pressure and 
shear stress relaxation in 2D is of order a factor of ten; in 
3D it is closer to a factor of four. Note that at higher tem- 
peratures (0.5-0.6) the relaxation time is of the same order 
as the discretization interval At — 20, so a large part of the 
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FIG. 7: Arrhenius plot of the inverse short-time rate — At/Ko 
determined from different fits to F{qm,t), compared with the re- 
laxation time Ta, and the cut-off time Ts of the 4-parameter fit, for 
2D (left panel) and 3D (right panel) systems. The data labelled "3- 
parameter (cons.)" refers to simultaneous fits of the three correlation 
functions at a given temperature, constraining the exponent a to be 
the same in all three. 
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FIG. 8: (Color online) Comparison of inverse short- time rate Ts = 
At/Ko for different kinds of correlation functions. Left and right 
panels show 2D and 3D data, respectively. Upper panels show the 
results of fitting each curve separately; lower panels show the results 
of constraining the power-law exponent a to be the same for all three 
curves at a given temperature. 



relaxation actually takes place within the first interval. Un- 
surprisingly, this limits the ability of the fitting procedure to 
accurately determine the zero-time value which leads to errors 
in normalization and hence in Ta, particularly for the collec- 
tive correlation functions of pressure and shear stress where 
the noise is high. In these cases a fit to a stretched exponential 
was made first, which was then used for the memory function 
fit. We note finally that is apparently nothing special happen- 
ing around the mode-coupling temperature Tc = 0.435 previ- 
ously identified for the 3D Kob-Andersen system (albeit with 
a slightly different cut-off in the potential). 



B. Short-time rate Ko and corresponding time Ts 
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We now investigate the parameters determined by the fit- 
ting procedure starting with the temperature dependence of 
the parameter Kq, related to the short-time relaxation rate. 
More precisely, we consider Ts = At/Ko, which is an ef- 
fective time scale quantifying the amount of relaxation that 
occurs over the interval At, a "short-time relaxation time". 
This is plotted Figure [T] along with for , t). Results 
of three different fitting schemes are shown: three-parameter 
fit, four-parameter fit (i.e., including the exponential cut-off) 
and a three-parameter fit where the exponent a is constrained 
to be equal for the three correlation functions at a given tem- 
perature. The motivation for the latter will become clear later 
on. The variation between the different fits gives an estimation 
of the error bars on this quantity, is less than r^, which is 
necessarily the case if < for m > 0. The main point is 
that Ta is more or less equal to the Tg at high temperatures — 
something necessarily true for exponential relaxation — ^but in- 
creases relative to it as temperature decreases. The tempera- 
ture dependence of the Tg is activated, which is to say it is at 



FIG. 9: (Color online) Temperature dependence of the exponent a 
for different quantities in 2D (left panels) and 3D (right panels) sim- 
ulations. Upper panels show the result of independently fitting the 
different correlation functions; lower panels show the common ex- 
ponent obtained by constraining it to be the same for all three, and 
that obtained in the (unconstrained) four-parameter fits to F[qm,t). 
In both systems, T =0.6 values are from fits using At =10. The 
horizontal lines indicate a =1.5. 



least Arrhenius. The data in the figure are not precise enough 
to determine whether the rate is super- Arrhenius, but certainly 
the effective barrier (the slope in the figure) is lower than for 
Ta itself. Figure [8] shows comparison of Tg for the three dif- 
ferent quantities. The upper panels represent independent fits, 
while the lower ones represent constrained fits where a is the 
same for the three quantities. 
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C. Power-law exponent a 



Figure |9] shows the temperature dependence of the power- 
law exponent a for the two different systems, the three dif- 
ferent correlation functions, and for different fitting schemes. 
Considering first the top two panels, the interesting feature 
is that the exponent values seem to converge as the temper- 
ature is lowered, towards ^1.6 in 2D and ^1.5 in 3D. At 
higher temperatures there is significant scatter, and, for the 
3D case at least, generally higher values up to around 2. The 
apparent convergence at low temperatures suggests that per- 
haps the true values of the exponents are in fact equal for the 
three functions at each temperature; the scatter in the uncon- 
strained fits would then be understood to be due to noise in the 
data. We therefore attempted a fit where a is constrained to 
be the same for all three correlation functions. In most cases 
the quality of the fit is indistinguishable by eye from the un- 
constrained fit (not shown). The exceptions are some fits to 
F{qm, t) in 2D where the constrained fit tends to reduce the 
value of a compared to the unconstrained fit, which results in 
a somewhat slower decay at long times (this could be proba- 
bly compensated for by including the exponential cut-off, as 
explained in the following paragraph; we have not done this). 
The temperature dependence of the constrained-a is more or 
less similar to the constrained case for the 3D data, while 
noticeably reduced in the 2D case. The effect on the con- 
strained fit on the values of can be seen in the lower panels 
of Fig. [8] There is a tendency towards smoother temperature 
dependence, which supports the idea that a constrained fit can 
yield more accurate results because the tendency to fit noise in 
any one of the curves is limited by the constraint. On the other 
hand the temperature dependence of Ts for the intermediate 
scattering function is actually less smooth in the constrained 
fit it has been "infected" by the greater noise in the collective 
functions. In principle this could be compensated for assign- 
ing a greater weight to the intermediate scattering function due 
to the smaller noise. This has not been attempted; it would 
require an unbiased estimate of the errors on the correlation 
functions. 



Also shown in the lower panels of Fig. |9] are values of a for 
unconstrained 4-parameter fits to F{q„i, t). What is notewor- 
thy here is that while values at higher temperatures still show 
appreciable scatter, there seems to be faster convergence at 
low temperatures, to values close to 1.55 and 1.50 in 2D and 
3D respectively (see in particular the four lowest temperature 
points in the 3D case). This suggests the interesting possibil- 
ity that the exponent could be in fact be independent of tem- 
perature as well as of which function is considered, and inde- 
pendent of which of the two systems. The apparently higher 
values at higher temperatures could be due to not including 
the exponential cut-off, which, like a higher value of a, would 
induce faster relaxation than a pure a = 1.5 power law. An 
attempt to find a temperature-independent a will be described 
below; first we consider the results for the parameter /. 
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FIG. 10: (Color online) Temperature dependence of the parameter / 
for different quantities in 2D (left panels) and 3D (right panels) simu- 
lations. Upper panels show results of independently fitting individual 
correlation functions; lower panels show the results of constraining 
a to be equal for all three. 



D. Tail-weight parameter / 

FigurefTOlshows the temperature dependence of the param- 
eter /, the weight of the negative tail in the memory func- 
tion relative to Kq (with no cut-off included). As T de- 
creases, / increases from a low value, approaching a value 
close to unity at the lowest temperatures. This increase ac- 
counts for that of the overall relaxation times though the 
factor 1/(1 — /). The value unity cannot be reached while 
remaining in equilibrium — because this implies a diverging 
relaxation time — unless the power-law behavior is cut off at 
the longest times. Thus there are two possibilities for the low- 
temperature behavior of /. The first is that it bends over and 
never reaches unity at finite temperature. In this case we can 
avoid including a cut-off, but the exact temperature depen- 
dence of / determines that of r^, given that a seems to be- 
come temperature-independent at such temperatures. The sec- 
ond possibility is that / becomes unity at a finite temperature, 
for example T ^0.3 in 2D or T ^0.4 in 3D. If this happens a 
cut-off must be included to keep finite. When considering 
the full temperature range there does seem to be a bend over 
towards smaller slopes at lower temperatures, but the data in 
Fig. [To] are not clean enough to draw a firm conclusion about 
the limiting T-dependence of /3 



E. Comparison with KWW fits, identification of common a 

It is interesting to compare with a more common charac- 
terization of the shape, the stretching parameter /3kww in 
the KWW function. As explained in Sec. IVIBI this may be 
compared to 2 — a, shown in Fig.[TT] The values of Pkww 
decrease as T decreases, while 2 — a increases. In a broad 
sense it seems that all are converging to a value around 0.5, 
but the degree of convergence for (3kww is much weaker. 
Considering the 3D data in particular, the values of 2 — a are 
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FIG. 11: (Color online) Comparison of 2 — a (determined using 3- 
parameter, unconstrained fits) with the stretching parameter /3kww 
in a fit to the stretched exponential for 2D (left) and 3D (right) sim- 
ulations.. Open (black) symbols are values of Pkww, taking values 
above 0.5, solids (red) symbols are values of 2 — a, taking values 
below 0.5. 
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FIG. 12: (Color online) Results of a constrained 4-parameter fit of 
F{qm,t) for all temperatures, where a is held fixed. Left panel 
shows data for 2D, right panel for 3D. 



within the range 0.47-0.52, while those of (3kww are spread 
over the range 0.55-0.7. The same trend can be seen, although 
less clearly, in the 2D data. It is possible that the Pkww 
values will converge at longer time scales, but the fact that 
a (or 2 — a) seems to converge more rapidly to a common 
value suggests that this representation may be more physi- 
cally meaningful. The value Pkww =0.5, corresponding to 
a =1.5, is interesting, because it has been argued theoretically 
and experimentally that this value is generic for structural "al- 
pha" relaxation in viscous liquids, see Refs. I20II2TII22II23II24I 
and references therein. 

Given the apparent convergence to a common value of a, 
and the hint in the 4-parameter fits to the ISF in 3D (lower- 
right panel of Fig. |9|l that this parameter may in fact be tem- 
perature dependent, we have tried to fit the high quality ISF 
data with a common value of a over the whole temperature 




FIG. 13: (Color online)Temperature dependence of Ts (symbols + 
solid lines) and Tc (symbols + dashed lines) for all correlation func- 
tions where a has been fixed at 1.58 (2D, left panels) and 1.5 (3D, 
right panels). While the 2D data does not fix Tc very well (probably 
because simply it is larger relative to the relaxation time and there- 
fore much more subject to noise in the tail), it seems plausible that a 
common value of Tc for the three correlation functions could apply 
in 3D. 



range, with the exponential cut-off included. The results are 
shown in Fig. [12] The common values of a that emerge from 
the fit are 1.58 for the 2D system and 1.50 for the 3D system. 
The quality of the "temperature-constrained" fits is is virtu- 
ally indistinguishable by eye from that of the unconstrained 
4-parameter fits (not shown). The latter, are of course, numeri- 
cally superior (since they are unconstrained), but this seems to 
be mostly due to better matching of the noise in the tail. This 
result strengthens the case for a having a true value of 1 .5 in 
3D, independent of temperature. The other interesting feature 
of this constrained-a fit is that the temperature dependence of 
the time-scale Tc of the exponential is much smoother, appar- 
ently tracking that of the relaxation time. What about the shear 
and pressure data? The data was of too low quality for an un- 
constrained four-parameter fit to work, but by fixing a to be 
equal to 1.5 we can include the exponential cut-off, making it 
effectively a three-parameter fit. The goodness of fit is consis- 
tently better than the original three-parameter fit (not shown), 
and now much smoother values, as well as Tc values are 
available for shear and pressure. These are shown in Fig.[T3] 
This works better for the 3D data than for the 2D data, which 
we suggest is due to the larger difference between and Tc 
in the latter case; the exponential cut-off has a significant ef- 
fect only at longer times where the signal to noise ratio in the 
coiTelation function is smaller, and the fit thus gives eiTatic 
values. 

In 3D the ratio Tc/tq, is in the range 2-3; in 2D the range 
is 8-9. Whether the cut-off time actually determines r^, or 
vice versa, or whether they determine each other in a self- 
consistent way, cannot be answered without a theoretical un- 
derstanding and/or an explicit model for the relaxation. Nor 
do we know what may be the cause of the different ratios in 
2D and 3D. To get a clearer picture more systems should be 
analyzed, and the possibility that Tc is size dependent cannot 
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FIG. 14: (Color online) Shear stress correlation functions and fits for 
2D system at T = 0.34. Main figure, comparison of true, inherent 
and fitted functions. The inherent function was based on minimiza- 
tion at intervals of 10 (1000 time steps); it has been rescaled to match 
the full correlation at t = 20 (by a factor 0.97). The inset shows a 
zoom in of the inherent correlation function and the 4-parameter fit 
with a = 1.58 at short times. In this case the inherent correlation 
function was based on minimization every 1 unit (100 time steps); 
here the rescale factor is 0.99. The need to rescale simply reflects the 
sampling error in the variance (t — value); that of the autocovari- 
ance is almost the same (for relatively short lags), hence different 
sampling results in curves which are more or less proportional to 
each other. When comparing normalized correlation at, say t = 20, 
the difference between minimization every 10 time steps (relatively 
short simulation times) and every 1000 time steps (relatively long 
simulation times) is of order 1%. 



be excluded. 



F. Comparison of the fitted Ko correspond witli tlie actual 
sliort-time decay of tiie inlierent quantities 

A feature of our fitting procedure is that we leave the zero- 
time value of the correlation function as a fitting parameter, 
since the actual zero-time value of the true correlation func- 
tion includes contributions from vibrational motion which 
rapidly decays. It is worth asking whether the value returned 
by the fit coincides with the value obtained by a careful deter- 
mination of the inherent correlation function. Since we focus 
on time intervals of order At ~ 20, we can run relatively long 
simulations, and obtain good statistics for inherent quantities, 
by minimizing at similar intervals. For example minimizing 
once every 1000 time steps (100 times less frequently than be- 
fore) corresponds to an interval At = 10, and means that the 
simulation time is not dominated as much by the cost of mini- 
mizing. This was done in both the 2D and 3D systems for two 
temperatures: 0.34 and 0.40 for 2D and 0.44 and 0.52 for 3D. 
In addition runs with minimization every 100 time steps were 
carried out for 2D, T = 0.34. 

TableUshows values of the normalized correlation function 
after one time interval At — 20, for two selected tempera- 
tures in each system, determined directly from the inherent 
state dynamics and inferred from different fits to the data for 



system 


T 


corr. 






V'(At)3-p 


V'(At)sE 


2D 


0.34 


ISF 


0.86 


0.94 


0.93 


0.95 


2D 


0.34 


shear 


0.64 


0.83 


0.74 


0.84 


2D 


0.34 


press. 


0.90 


0.95 


0.95 


0.95 


2D 


0.40 


ISF 


0.72 


0.83 


0.71 


0.84 


2D 


0.40 


shear 


0.46 


0.56 


0.60 


0.54 


2D 


0.40 


press 


0.77 


0.85 


0.70 


0.82 


3D 


0.44 


ISF 


0.94 


0.97 


0.97 


0.98 


3D 


0.44 


shear 


0.78 


0.91 


0.90 


0.89 


3D 


0.44 


press. 


0.91 


0.94 


0.93 


0.96 


3D 


0.52 


ISF 


0.70 


0.79 


0.73 


0.81 


3D 


0.52 


shear 


0.34 


0.43 


0.46 


0.40 


3D 


0.52 


press. 


0.65 


0.68 


0.66 


0.73 



TABLE I: Value of correlation functions (normalized to equal unity 
at time zero) at time At — 20; comparison between actual inherent 
value ("I") determined from the autocorrelation function of inherent 
quantities and values inferred from fitting the correlation function at 
longer times. Values for three fits are shown: ("1.5") the fit with a 
fixed at be 1.5, including the exponential cut-off; ("3 — p") the 3- 
parameter fit with no cut-off; and ("SE") the stretched exponential 
fit. The correlation at At implied by the fits is higher than the actual 
correlation; alternatively put, the amount of decorrelation that takes 
place on the time scale At is less than would be expected from the 
fits. 



t > At. Although there are some exceptions (particularly 
for 2D, T = 0.40), the trend is that the value of the normal- 
ized inherent correlation implied by the fits is higher than the 
actual value. For example for 2D at T = 0.34, considering 
the intermediate scattering function, the actual correlation is 
0.86, while the fits all imply a value between 0.93-0.94. An- 
other way of putting this is that the amount of decorrelation 
that takes place between time zero and time At is greater than 
is implied by the fits of the correlation at longer times (note 
that tpi'^t) = "01 1 ^ Kq). This can be thought of as 
somewhat analogous to the extra relaxation associated with 
the vibrational motion which has been removed, but this is 
part of the inherent relaxation. Schr0der et al. also noticed 
a short-time component of the inherent structural relaxation 
(intermediate scattering function) which was not accounted 
for by stretched exponential fits.'^ This was particularly no- 
ticeable at high temperatures and appeared to disappear in the 
activated regime. The comparison above suggests that this ex- 
tra relaxation does not vanish, even at quite low temperatures. 
This discrepancy, between the very short-time relaxation and 
that for later times (t > At), presumably corresponds to the 
change of apparent power-law exponent in the memory func- 
tion around this time scale (section HVBI ). One could specu- 
late that the existence of this extra component of the inherent 
relaxation could be associated with very small energy barriers, 
much smaller than the temperature. Such barriers would sepa- 
rate distinct inherent states, and thus "transitions" would show 
up in the inherent correlation functions, but these transitions 
would not be activated. The author has determined that the 
energy barriers in these systems are exponentially distributed, 
thus there are always some small barriers*^ 
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FIG. 15: (Color online) Wavenumber-dependent intermediate scat- 
tering function (symbols and dotted lines) and 4-parameter fits with 
a = 1.5 (solid lines), for 3D system, r=0.48. The inset shows the 
g-dependence of the tail-weight parameter /. 
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FIG. 16: (Color online) Main figure, double-log plot of the 
g— dependence of the time-scales associated with the 4-parameter fit: 
the relaxation time r^, the inverse of short-time rate, Ts, and the ex- 
ponential cut-off time Tc. A line proportional to is shown as a 
guide to the eye. The q values used are as in Fig. 1151 



G. Wavenumber dependence of intermediate scattering 
function 

For simplicity, in studying the intermediate scattering 
above, we have restricted the analysis to a particular 
wavenumber, q™, and to large particles. Given that a common 
value of the exponent a seems to be valid for the three correla- 
tion functions we have focused on so far, it makes sense to ask 
whether this applies for the intermediate scattering function at 
other g-values. In particular, what happens at low q, where we 
know that the relaxation becomes more exponential? Fig.fTS] 
shows data and fits for 3D at T = 0.46. Here 4-parameter fits 
with fixed a = 1.5 were used. 

Figure [16] shows the ^-dependence of the relaxation time 
Tq, inverse of the short time rate, Ts and the cut-off time of the 
memory function t^. The inset shows the (/-dependence of /. 
While the power-law exponent is fixed at 1.5 for all q-values, 
/ varies strongly, approaching zero at small q. This corre- 



sponds to relaxation that becomes exponential in this limit. 
Consistent with this, and Tq approach each other Tc seems 
to follow Tot- In normal diffusion one expects the relaxation 
time to be proportional to g^^, indicated by the dashed line. 
The relaxation time seems to approach this dependence in the 
low-(7 limit (as is expected, given that diffusion should be nor- 
mal at long length and time scales). More interesting is that 
Ts seems to show an inverse square dependence over most of 
its range, corresponding to normal-diffusive behavior being 
observed at short times. This will be discussed below. 



VI. DISCUSSION 

A. Generalized Langevin equations and continuous-time 
random walks 

What does a negative inverse power-law memory function 
mean? One interpretation derives from the formalism of the 
generalized Langevin equation (GLE), which offers a frame- 
work for interpreting memory functions. The GLE for a dy- 
namical variable A{t) is a stochastic differential equation 

dA /■* 

— = K{t-t')A{t')dt' + Fit). (17) 
d-t Jq 

The first term on the right hand side involves an effective 
friction, non-local in time, where the memory kernel K{t~t') 
connects the friction at current time t with the value of A at 
a (previous) time t' . The last term F{t) is the so-called ran- 
dom force, or noise term. According to the second fluctuation- 
dissipation theorem, the memory kernel is the autocorrelation 
function of the noise, or 

K{t) = {F{0)F{t)). (18) 

It may also be shown that this memory kernel K{t) is 
the memory function associated with the autocorrelation of 
A{t) — the same memory function used in the analysis of this 
paper Therefore one can interpret the obtained memory func- 
tions as the autocorrelation functions of the random force of 
the corresponding GLE. By "force" is meant the stochastic 
increment to the quantity A; in the original version of the 
Langevin equation A was the velocity of a particle, whose 
change was indeed a force (divided by mass). Negative K{t) 
for t > implies that the changes, which are presumably re- 
lated to the so-called "flow-events" which define the long-time 
dynamics, are anti-correlated. The simplest interpretation of 
this would be forward-backward correlation of events — that 
is, an event is more likely to be the reverse of the previous one 
than to be something else. The parameter / is then a direct 
measure of this anti-correlation. As it approaches unity, this 
tendency becomes so strong that the relaxation time diverges 
(see the next subsection). 

An alternative interpretation involves uncorrected events 
but with a non-exponential waiting-time distribution. This can 
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be understood within the continuous-time random walk for- 
malism of Montroll and Weiss. The Montroll-Weiss equa- 
tion relates the Laplace transform of F{q,t) to the Laplace 
and Fourier transforms of the waiting-time and jump-size dis- 
tributions, p^(s) and pj{k), respectively. It is straightforward 
to use it to derive a relation between the waiting-time distri- 
bution and the memory function. Skipping the details, this 
yields 

p^{s) = ^ . (19) 

This relation can be considered a formal equivalence of 
the two pictures — anti-correlated random forces and non- 
exponential waiting-time distribution. The question of which 
picture is physically more relevant can only be determined by 
more detailed analysis of the dynamics in simulations. The 
existence of forward-backward correlations has indeed been 
noted by Heuer and coworkers (see Ref. 127" for a review of 
this work), but they argue that these are trivial, operating only 
on relatively short time scales. They have proposed a proce- 
dure for coarse-graining the dynamics by grouping the basins 
of attraction of inherent states together to form "meta-basins" 
(MBs). The definition of an MB is based on the requirement 
that once the system has escaped from an MB, it is unlikely 
to return there. They have shown that the dynamics on long 
time scales is essentially a random walk among MBs, but with 
non-trivial waiting-time distribution.^^ If we accept this inter- 
pretation, the memory function analysis presented here could 
be used as a way to find the waiting-time distribution from the 
correlation function. For example, in the / ^ 1 limit with 
no cut-off (see next section) one finds that the waiting time 
distribution has a power-law tail with the exponent a as the 
memory function itself. 

On the other hand, is it possible that the equivalence be- 
tween the waiting time distribution point of view and the 
forward-backward point of view is in fact more than formal? 
The non-exponential distribution of waiting times could seen 
as a consequence of forward-backward correlations. The size 
of an MB is not a fixed property of the landscape^^ instead be- 
ing defined in terms of dynamical behavior: a given MB can 
increase in size as temperature decreases, reflecting the in- 
creased tendency for returns (the system needs to get further 
away on average before its chance of escaping for good ex- 
ceeds 50%). Correspondingly, the "anti-correlation tendency" 
/ increases towards unity; as long as it is less than unity, how- 
ever, a diffusive regime is reached — the system will eventu- 
ally forget where it started. In the "CTRW picture", as long as 
the waiting-time distribution has a finite mean, upon sufficient 
coarse-graining it will become exponential and therefore the 
dynamics will be diffusive in the normal sense. The point of 
Heuer and co-workers is that there is a coarse-graining time 
scale at which events are uncorrected but the waiting-time 
distribution is non-trivial (non-exponential). It is also true 
in the "GLE picture" involving anti-correlated events that, 
as long as / < 1 or there is an exponential cut-off, a nor- 
mal diffusive regime must also eventually be reached under 
coarse-graining in time. But it is not clear yet whether coarse- 



graining in the GLE picture can produce a dynamics that cor- 
responds to the CTRW picture. This is a matter for future 
investigations. 



B. The f — 1 limit: sub-diffusion 

In this section we discuss briefly the limit / ^ 1, in the 
case where there is no exponential cut-off. As mentioned, in 
this limit the relaxation time diverges, but the dynamics de- 
fined by this memory function are still meaningful, and have 
been studied in the context of anomalous diffusion, fractional 
Brownian motion and related processes. Applications in- 
clude protein dynamics-'^'^-' and dielectric response of glassy 
medial^ A noise term in a GLE with this kind of correlation 
function is termed fractional Gaussian noise; when it acts on 
an unbound degree of freedom, the latter undergoes fractional 
Brownian motion. Correlation functions generally involve a 
function known as the Mittag-Leffler (ML) function:^^ They 
behave at intermediate times like a stretched exponential, with 
stretching parameter (3kww = 2 — a, before switching to a 
power-law at long times, such that the integral is infinite. This 
is why it made sense to compare 2 — a with the fitted (3kww 
in Fig. [m 

There are two broad classes corresponding to positively or 
negatively correlated noise. When comparing to the litera- 
ture on fractional diffusion and related processes it is impor- 
tant to realize that the Langevin equation is typically formu- 
lated with the dynamical variable A{t) as the velocity of a 
particle. For a free particle, if the noise is fractional Gaus- 
sian with negative correlation, then the velocity autocorrela- 
tion function is an ML function. Its integral diverges, there- 
fore the position dynamics is actually super-diffusive. Appli- 
cations more typically involve sub-diffusive motion of a har- 
monically bound particle.''^ The noise acting on the velocity 
is positively correlated and it is the position autocorrelation 
function that decays as an ML function. 

The memory function corresponding to this limit is K{s) = 
s"~^, which is clearly scale invariant: Rescaling, or coarse- 
graining in time produces the same dynamics. This corre- 
sponds to a fixed point in the renormalization group sense. 
For negatively correlated noise, the fixed point is unstable 
(see the chapter by Qian in Ref. l30l) — ^if we are not quite in 
this limit, then rescaling will eventually takes us away from it 
and back to the normal diffusive regime, as discussed in the 
previous subsection. The results of the analysis presented in 
this paper could be therefore interpreted as saying that glassy 
dynamics approaches sub-diffusive dynamics with a particu- 
lar exponent a, which (in 3D) is very close to 1.5. As the 
limit / = 1 is approached the relaxation becomes close to 
a Mittag-Leffler function, resembling a stretched exponential 
with /3kww = 0.5. As we have seen, even at the lowest tem- 
peratures, actually fitting with a stretched exponential yields a 
different exponent, so it is hard to see the convergence to 0.5, 
while this is clearer in the memory function based approach. 
The theoretical challenge is, of course, to understand where 
this exponent comes from, although there are models for vis- 
cous liquid dynamics which predict it (Refs. I20ll2lll22i23ll24l 
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and references within). 

Does this limit ever get reached? Although / clearly ap- 
proaches unity, it has not reached it yet at the temperatures 
simulated so far Moreover, there is always an exponential 
cut-off in the memory function. Thus regardless of whether 
/ actually reaches unity or not, it will at some point be the 
exponential cut-off which will determine the relaxation time. 
Looking at Fig.[T3] for example, suggests that this will happen 
around T = 0.40 in 3D. Simple-minded linear or quadratic 
extrapolation of / suggests this also. If we set / = 1 for 
temperatures below say T — 0.40 and assume that the expo- 
nential cut-off, Tc is a fixed multiple of Tq, (actually the ratio 
seems to increase), then one can make an extrapolation of Tq, 
which gives a sharp crossover to an almost Arrhenius form 
(not shown). If / approaches unity more smoothly, then the 
crossover will not be so abrupt, but in any case we predict 
that a change in the temperature dependence will occur around 
T = 0.40. 



C. Short-time dynamics 

The parameter = At/K^ is a measure of the dynamics 
at short times. By fixing the time scale, it measures how much 
activity happens in a given time, as opposed to standard mea- 
sures of relaxation, which concern how much time is required 
for a given amount of relaxation to occur. For exponential re- 
laxation there is no difference, for relatively small intervals 
(or long relaxation time Tc): If the interval is At, then the 
amount of correlation at this time is 



1 - Ko = CXp{-At/Ta) ^1 - At/Ta, (20) 

which gives Tg — Tq. Thus a difference between these two 
quantities is equivalent to non-exponential relaxation. When 
At is much smaller than Tq,, the expected number of activated 
events that take place in a small volume in a time interval of 
this length is much smaller than unity, thus we cannot expect 
to see any correlation effects. We suggest therefore that is 
like a "bare" relaxation time, measuring the dynamics before 
correlation effects play a role. Recall (Fig.[T4]i that there is a 
noticeable amount of (inherent) relaxation in this time inter- 
val which is not accounted for by the fitted Tg. We assume this 
to be non-activated relaxation, and not relevant for what takes 
place at longer time scales (in this sense it resembles the vibra- 
tional part of the relaxation). The evidence supporting the idea 
of Tg as a measure of the bare, uncorrected relaxation rate is 
the near- Arrhenius temperature dependence (see Fig.[T3]l, and 
the near wavenumber dependence (Fig. [TST i The use of 
a fixed, relatively short time interval to extract an apparent 
"underlying Arrhenius behavior" in viscous liquid dynamics 
was recently demonstrated by de Souza and Wales^^'^^. They 
studied the diffusion constant determined by the mean squared 
displacement in a fixed time interval of 25 LJ units. They 
were also able to explicitly relate the deviations from expo- 
nential relaxation to anti-correlation of events in subsequent 
time windows. The present results are essentially equivalent. 



but obtained in the more general context of studying an arbi- 
trary autocorrelation function. 



VII. CONCLUSION 

We have presented a method for fitting the kind of slowly 
decaying autocorrelation functions typical of the dynamics of 
glass-forming liquids. An explicit functional form is postu- 
lated, not for the autocorrelation function, but for its asso- 
ciated memory function: a positive zero-time parameter Kq, 
and a negative inverse power law for non-zero time. An long- 
time exponential cut-off may optionally be included. This in- 
volves another parameter, although we have shown that good 
fits may be obtained by fixing the exponent to a common 
value, reducing the number of parameters per data set essen- 
tially to three. While it is possible to obtain the exact memory 
function from the autocorrelation function, it is not useful to 
fit the former directly. Rather we match the autocorrelation 
function by optimizing the parameters defining the memory 
function. Using a quite general functional form for K{t) we 
noted a crossover from one apparent power-law to another at 
relatively short times, which allowed a simplification by con- 
sidering data separated at times longer than the cross-over, 
namely the simple power-law form (with possible exponential 
cut-off). In particular we chose At = 20 in Lennard-Jones 
units for all of the analysis. While the fits are about as good as 
fits with a stretched exponential, the parameters obtained have 
arguably greater physical significance than those of the latter. 
In particular the exponent a seems to be more or less inde- 
pendent of the relaxing quantity, and indeed of temperature 
(as long as the exponential cut-off is included), which is not 
the case for the stretching parameter Pkww- Moreover the 
parameter Ko may be interpreted as a short time rate; its tem- 
perature dependence is, interesting, significantly weaker than 
that of the alpha relaxation times for the different autocorrela- 
tion functions. As temperature decreases, the amplitude of the 
power tends to approach a limiting value associated with the 
mathematical description of anomalous diffusion/relaxation 
characterized mathematically by the Mittag-Leffler function. 
We hypothesize that this value is reached at finite temperature 
but the associated divergence of relaxation time is removed 
by the presence of an exponential cut-off. Including the lat- 
ter in the fits for the self-intermediate scattering function al- 
lowed the identification of a third time scale, longer than the 
alpha time. It is hypothesized that the strong non-Arrhenius 
temperature dependence is due to a crossover from the time 
associated with the short-time rate to that associated with the 
long-time cut-off of the memory function, and that in particu- 
lar at temperatures somewhat lower than those simulated, the 
non-Arrhenius behavior will weaken noticeably. 
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APPENDIX A: MEMORY FUNCTION FITTING 

Here we describe in a little more detail how we optimize 
the fit. We start by describing the procedure using the general 
fitting form which includes the piecewise linear function L{t). 
Suppose first we have a correlation function based on inherent 
dynamics, so that there is no vibrational decay at short times, 
and in particular the t = value is known. Then we may 
normalize by the latter to get the normalized autocorrelation 
function ip{t) — ipn where n — t/ (At) n is the integer corre- 
sponding to t. We write the fitting form as Kl'^*{{pa}), where 
pi,P2, ■ ■ ■ are the parameters of the fitting form. The node 
values Lk are not used directly as parameters; rather their log- 
arithms are. This constrains them to be positive and minimizes 
numerical problems associated with the broad range of values 
that they turn out to have. For given values of the {pa}, we can 
easily compute the associated correlation function {{Pa}) 
and then the objective function 

FobjiiPa}) = ~ 'f'senit, {Pa})f ■ (Al) 

n 

To minimize Fobj we used the conjugate gradient techniquer^ 
calculating the gradient numerically. 

For the main data analysis we switched to using simpler 
functional forms without the piecewise linear function. Also 
significant is that we used only the true dynamical correla- 
tion functions (as opposed to those of the inherent dynamics) 
which are more accurate and smooth; this removed the need 
to do expensive minimizations. By considering times at inter- 
vals of order 20, we could avoid the vibrational decay, except 
that now the correct initial value of correlation function had 
to be treated as an unknown parameter Since this affects the 
normalization, it is important to handle this carefully. Sup- 
pose we have a non-normalized correlation function at dis- 



crete times, not including time zero: Ci, C2, C3, . . .. Since 
we are going to normalize anyway, we choose as the param- 
eter not Co, but the ratio Rqi = Cq/Ci. And since this 
is required be greater than unity, it is represented internally 
(i.e., within the fitting algorithm) in terms of a logarithm Or: 
Roi = 1 + exp(^/{), where Or is the actual parameter varied 
by the algorithm. Note that Rqi is directly related to Kq: 

ifo = l-V'i = l-(Ci/Co) = l-l/i?, (A2) 

so that one may equivalently consider Kq as the fitting param- 
eter An important technicality must be mentioned here. If the 
objective function is defined as above, there is a problem due 
to the fact that the normalization factor is now an adjustable 
parameter. This means that the objective function could be re- 
duced by making Rqi, and hence Cq, larger. This was found to 
be a problem for data sets from higher temperatures where the 
relaxation times were relatively short and therefore the data 
was not sufficient to constrain Rqi to reasonable values. The 
problem was solved by normalizing by the fixed quantity Ci 
instead (but only for the purposes of defining the objective 
function). For n > (times greater than A< = 20) we used 
in some cases the single power law of Eq. [14] The tail-weight 
fraction / is constrained to lie between and 1 by expressing 
it internally as the / = arctan(/i) where /i is unconstrained. 
The exponent a was represented as itself because it tended to 
vary within a relatively small range (1.5-2.2), with no risk of 
assuming values which would cause numerical problems. In 
other cases an exponential cut-off factor was included, Eq.fTSl 
which involves the timescale Tcut- Because this can vary over 
several orders of magnitude, to aid convergence it was ex- 
pressed internally in terms of its logarithm, Tcut ~ exp{109r), 
where the factor 10 was found to improve convergence (pre- 
sumably by making the typical values and range of 9r compa- 
rable to those of other parameters). 
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We can exclude Arrhenius dependence (of 1 — /), as this would 
imply Arrhenius dependence of Tq, 

The argument of the Mittag-Leffler function is actually (t/rA/i, ) 
for a time scale tml; so the correlation function should really 
be called a "stretched Mittag-Leffler function", but it seems to 
be standard to refer to the correlation function also as a Mittag- 
Leffler function, which is not strictly correct. 
These can also be defined in terms of the so-called Hurst exponent 
H being respectively greater than i, which also corresponds to 
< a < 1 or less than i, corresponding to 1 < a < 2. At 
iif = i the exponent is formaUy equal to unity but its coefficient 
vanishes thus white (uncorrelated) noise is recovered. 



